
<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta charset="utf-8" />
    <title>psf &#8212; Texas A&amp;M Oil spill / Outfall Calculator 2.1.0 documentation</title>
    <link rel="stylesheet" href="../_static/alabaster.css" type="text/css" />
    <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
    <script id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
    <script src="../_static/jquery.js"></script>
    <script src="../_static/underscore.js"></script>
    <script src="../_static/doctools.js"></script>
    <script src="../_static/language_data.js"></script>
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
   
  <link rel="stylesheet" href="../_static/custom.css" type="text/css" />
  
  
  <meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />

  </head><body>
  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          

          <div class="body" role="main">
            
  <h1>Source code for psf</h1><div class="highlight"><pre>
<span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Particle Size Functions</span>
<span class="sd">=======================</span>

<span class="sd">This module provides empirical functions for computing particle sizes (bubble</span>
<span class="sd">and droplets) from jet releases into water.  Model equations are from::</span>

<span class="sd">    Johansen et al. (2013)</span>
<span class="sd">    Li et al. (2017)</span>
<span class="sd">    Wang et al. (2018)</span>

<span class="sd">Utilities are also included to generate log-normal or Rosin-Rammler size</span>
<span class="sd">distributions and to apply the d_95 rule (the idea that the d_95 value of a</span>
<span class="sd">volume-size distribution should not exceed the maximum stable particle size).</span>

<span class="sd">The functions in this module are used by the `particle_size_models` module</span>
<span class="sd">and are not intended to be called directly.  In the future, these functions</span>
<span class="sd">could, for example, be ported to Fortran or another language and wrapped </span>
<span class="sd">in Python for use by the `particle_size_models` class objects.  </span>

<span class="sd">See Also</span>
<span class="sd">--------</span>
<span class="sd">particle_size_models.ModelBase, particle_size_models.PureJet, </span>
<span class="sd">particle_size_models.Model</span>

<span class="sd">Notes</span>
<span class="sd">-----</span>
<span class="sd">Particle size distributions are computed from the following sources::</span>

<span class="sd">* Johansen, O., Brandvik, P. J., and Farooq, U. (2013) &quot;Droplet breakup in </span>
<span class="sd">  subsea oil releases - Part 2: Predictions of droplet size distributions </span>
<span class="sd">  with and without injection of chemical dispersants.&quot; Mar Pollut Bull, </span>
<span class="sd">  73(1), 327-335. This reference is intended for liquid oil breakup only.</span>

<span class="sd">* Li, Z., Spaulding, M., French McCay, D., Crowley, D., and Payne, J. R. </span>
<span class="sd">  (2017) &quot;Development of a unified oil droplet size distribution model with </span>
<span class="sd">  application to surface breaking waves and subsea blowout releases </span>
<span class="sd">  considering dispersant effects.&quot; Mar Pollut Bull, 114(1), 247-257.  The</span>
<span class="sd">  authors apply this reference for breakup of gas bubble or oil droplets.</span>

<span class="sd">* Wang, B., Socolofsky, S. A., Lai, C. C. K., Adams, E. E., and Boufadel, M. </span>
<span class="sd">  C. (2018). &quot;Behavior and dynamics of bubble breakup in gas pipeline leaks </span>
<span class="sd">  and accidental subsea oil well blowouts.&quot; Mar Pollut Bull, 131, 72-86.  </span>
<span class="sd">  This reference is only intended for gas bubble breakup.</span>

<span class="sd">The maximum stable droplet size of an immiscible liquid in seawater is taken</span>
<span class="sd">from::</span>

<span class="sd">* Clift, R., Grace, J., and Weber, M. E. (1978) *Bubbles, Drops, and \ </span>
<span class="sd">  Particles*, Dover Publications, Inc., Mineola, New York.</span>

<span class="sd">For gas bubbles in water, the maximum stable bubble size is taken from a </span>
<span class="sd">method in Grace et al::</span>

<span class="sd">* Grace, J.R., Wairegi, T., Brophy, J., (1978) &quot;Break-up of drops and bubbles </span>
<span class="sd">  in stagnant media,&quot; Can. J. Chem. Eng. 56 (1), 3-8.</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="c1"># S. Socolofsky, December 2019, Texas A&amp;M University &lt;socolofs@tamu.edu&gt;</span>

<span class="kn">from</span> <span class="nn">__future__</span> <span class="kn">import</span> <span class="p">(</span><span class="n">absolute_import</span><span class="p">,</span> <span class="n">division</span><span class="p">,</span> <span class="n">print_function</span><span class="p">)</span>

<span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">dbm</span>

<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="kn">import</span> <span class="n">fsolve</span>
<span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="kn">import</span> <span class="n">minimize</span>

<span class="c1"># Use SI units throughout</span>
<span class="n">G</span> <span class="o">=</span> <span class="mf">9.81</span>

<span class="c1"># General Utilities ----------------------------------------------------------</span>

<div class="viewcode-block" id="mass2vol"><a class="viewcode-back" href="../autodoc/psd/psf.mass2vol.html#psf.mass2vol">[docs]</a><span class="k">def</span> <span class="nf">mass2vol</span><span class="p">(</span><span class="n">m</span><span class="p">,</span> <span class="n">rho</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert a mass or mass flux to a volume or volume flux</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    m : ndarray</span>
<span class="sd">        Array of masses (kg) or mass fluxes (kg/s) for each component of a </span>
<span class="sd">        given fluid</span>
<span class="sd">    rho : float</span>
<span class="sd">        In-situ density of a given fluid (kg/m^3)</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    q : float</span>
<span class="sd">        Corresponding volume (m^3) or volume flux (m^3/s) of a given fluid</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute volume handling zero-fluxes correctly</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">m</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">q</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">m</span><span class="p">)</span> <span class="o">/</span> <span class="n">rho</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">q</span> <span class="o">=</span> <span class="mf">0.</span>
    
    <span class="k">return</span> <span class="n">q</span></div>


<span class="c1"># Probabiity Density Functions -----------------------------------------------</span>

<div class="viewcode-block" id="rosin_rammler"><a class="viewcode-back" href="../autodoc/psd/psf.rosin_rammler.html#psf.rosin_rammler">[docs]</a><span class="k">def</span> <span class="nf">rosin_rammler</span><span class="p">(</span><span class="n">nbins</span><span class="p">,</span> <span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return the volume size distribution from the Rosin-Rammler distribution</span>
<span class="sd">    </span>
<span class="sd">    Returns the fluid particle diameters in the selected number of bins on</span>
<span class="sd">    a volume-fraction basis from the Rosin Rammler distribution with </span>
<span class="sd">    parameters d_50, k, and alpha.  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    nbins : int</span>
<span class="sd">        Desired number of size bins in the particle volume size distribution</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume mean particle diameter (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter of the Rosin-Rammler distribution (=log(0.5) for d_50)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter of the Rosin-Rammler distribution</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    de : ndarray</span>
<span class="sd">        Array of particle sizes at the center of each bin in the distribution</span>
<span class="sd">        (m)</span>
<span class="sd">    vf : ndarray</span>
<span class="sd">        Volume fraction in each bin (--)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the de/d50 ratio for the edges of each bin in the distribution </span>
    <span class="c1"># using a log-spacing</span>
    <span class="n">a99</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.995</span><span class="p">)</span> <span class="o">/</span> <span class="n">k</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">alpha</span><span class="p">)</span>
    <span class="n">a01</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.01</span><span class="p">)</span> <span class="o">/</span> <span class="n">k</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">alpha</span><span class="p">)</span>
    <span class="n">bin_edges</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">a01</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">a99</span><span class="p">),</span> <span class="n">nbins</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span>
    
    <span class="c1"># Find the logarithmic average location of the center of each bin</span>
    <span class="n">bin_centers</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">)</span>
    
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">)):</span>
        <span class="n">bin_centers</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">+</span> 
                         <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="p">]))</span> <span class="o">/</span> <span class="mf">2.</span><span class="p">)</span>
    
    <span class="c1"># Compute the actual diameters of each particle</span>
    <span class="n">de</span> <span class="o">=</span> <span class="n">d50</span> <span class="o">*</span> <span class="n">bin_centers</span>
    
    <span class="c1"># Get the volume fraction within each bin</span>
    <span class="k">if</span> <span class="n">d50</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">vn</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">k</span> <span class="o">*</span> <span class="n">bin_edges</span><span class="o">**</span><span class="n">alpha</span><span class="p">)</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">))</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">):</span>
            <span class="n">vf</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">vn</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">vn</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="n">vf</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">vf</span><span class="p">)</span>
    
    <span class="c1"># Return the particle size distribution</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">vf</span><span class="p">)</span></div>


<div class="viewcode-block" id="log_normal"><a class="viewcode-back" href="../autodoc/psd/psf.log_normal.html#psf.log_normal">[docs]</a><span class="k">def</span> <span class="nf">log_normal</span><span class="p">(</span><span class="n">nbins</span><span class="p">,</span> <span class="n">d50</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return the volume size distribution from the Log-normal distribution</span>
<span class="sd">    </span>
<span class="sd">    Returns the fluid particle diameters in the selected number of bins on</span>
<span class="sd">    a volume-fraction basis from the Log-normal distribution with parameters</span>
<span class="sd">    d50 and sigma_x.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    nbins : int</span>
<span class="sd">        Desired number of size bins in the particle volume size distribution</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume mean particle diameter (m)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    de : ndarray</span>
<span class="sd">        Array of particle sizes at the center of each bin in the distribution</span>
<span class="sd">        (m)</span>
<span class="sd">    vf : ndarray</span>
<span class="sd">        Volume fraction in each bin (--)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function uses the log-normal distribution defined at::</span>
<span class="sd">    </span>
<span class="sd">        https://en.wikipedia.org/wiki/Log-normal_distribution</span>
<span class="sd">    </span>
<span class="sd">    last accessed on 1/10/2020.  The relationship between mu_x and sigma_x </span>
<span class="sd">    of the real x values and mu and sigma of the log(x) values is from the</span>
<span class="sd">    notes for ENGR 102 posted here:</span>
<span class="sd">        </span>
<span class="sd">        https://ceprofs.civil.tamu.edu/ssocolofsky/ENGR102/</span>
<span class="sd">        Downloads/19c/Week_12/week_12.pdf</span>
<span class="sd">    </span>
<span class="sd">    last accessed on 1/10/2020.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the de/d50 ratio for the edges of each bin in the distribution </span>
    <span class="c1"># using a log-spacing</span>
    <span class="n">a0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d50</span><span class="p">)</span> <span class="o">-</span> <span class="mf">2.8</span> <span class="o">*</span> <span class="n">sigma</span><span class="p">)</span> <span class="o">/</span> <span class="n">d50</span>
    <span class="n">a1</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d50</span><span class="p">)</span> <span class="o">+</span> <span class="mf">2.3</span> <span class="o">*</span> <span class="n">sigma</span><span class="p">)</span> <span class="o">/</span> <span class="n">d50</span>
    <span class="n">bin_edges</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">a0</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">a1</span><span class="p">),</span> <span class="n">nbins</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span>
    
    <span class="c1"># Find the logarithmic average location of the center of each bin</span>
    <span class="n">bin_centers</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">)</span>
    
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">)):</span>
        <span class="n">bin_centers</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">+</span> 
                         <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_edges</span><span class="p">[</span><span class="n">i</span><span class="p">]))</span> <span class="o">/</span> <span class="mf">2.</span><span class="p">)</span>
    
    <span class="c1"># Compute the actual diameters of each particle</span>
    <span class="n">de</span> <span class="o">=</span> <span class="n">d50</span> <span class="o">*</span> <span class="n">bin_centers</span>
    
    <span class="c1"># Compute log-normal parameters for de/d50 distribution</span>
    <span class="n">mu</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span><span class="p">)</span>
    
    <span class="c1"># Get the volume fraction within each bin</span>
    <span class="k">if</span> <span class="n">d50</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="n">bin_centers</span> <span class="o">*</span> <span class="mf">1.</span> <span class="o">/</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">))</span> <span class="o">*</span> \
             <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">bin_centers</span><span class="p">)</span> <span class="o">-</span> <span class="n">mu</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">sigma</span><span class="o">**</span><span class="mi">2</span><span class="p">))</span>
        <span class="n">vf</span> <span class="o">=</span> <span class="n">vf</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">vf</span><span class="p">)</span>
    
    <span class="c1"># Return the particle size distribution</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">vf</span><span class="p">)</span></div>


<div class="viewcode-block" id="ln2rr"><a class="viewcode-back" href="../autodoc/psd/psf.ln2rr.html#psf.ln2rr">[docs]</a><span class="k">def</span> <span class="nf">ln2rr</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">sigma</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert the parameters of a log-normal distribution to Rosin-Rammler</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        The median particle size of a volume distribution</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume mean particle diameter (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter of the Rosin-Rammler distribution (=log(0.5) for d_50)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter of the Rosin-Rammler distribution</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute d95 of the log-normal distribution</span>
    <span class="n">mu</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d50</span><span class="p">)</span>
    <span class="n">mu_95</span> <span class="o">=</span> <span class="n">mu</span> <span class="o">+</span> <span class="mf">1.6449</span> <span class="o">*</span> <span class="n">sigma</span>
    <span class="n">d95</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">mu_95</span><span class="p">)</span>
    
    <span class="c1"># Find parameters of Rosin-Rammler with same d50 and d95</span>
    <span class="n">k</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">0.5</span><span class="p">)</span>
    <span class="n">alpha</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.95</span><span class="p">)</span> <span class="o">/</span> <span class="n">k</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d95</span> <span class="o">/</span> <span class="n">d50</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span></div>


<div class="viewcode-block" id="rr2ln"><a class="viewcode-back" href="../autodoc/psd/psf.rr2ln.html#psf.rr2ln">[docs]</a><span class="k">def</span> <span class="nf">rr2ln</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert the parameters of a Rosin-Rammler distribution to log-normal</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume mean particle diameter (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter of the Rosin-Rammler distribution (=log(0.5) for d_50)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter of the Rosin-Rammler distribution</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        The median particle size of a volume distribution</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute d95 of the Rosin-Rammler distribution</span>
    <span class="n">k</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">0.5</span><span class="p">)</span>
    <span class="n">d95</span> <span class="o">=</span> <span class="n">d50</span> <span class="o">*</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.95</span><span class="p">)</span> <span class="o">/</span> <span class="n">k</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">alpha</span><span class="p">)</span>
    
    <span class="c1"># Find the parameters of log-normal with same d50 and d95</span>
    <span class="n">sigma</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d95</span><span class="p">)</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d50</span><span class="p">))</span> <span class="o">/</span> <span class="mf">1.6449</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span></div>


<div class="viewcode-block" id="rosin_rammler_fit"><a class="viewcode-back" href="../autodoc/psd/psf.rosin_rammler_fit.html#psf.rosin_rammler_fit">[docs]</a><span class="k">def</span> <span class="nf">rosin_rammler_fit</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">d_max</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=</span><span class="mf">1.8</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return d_50, k, and alpha for the Rosin-Rammler distribution</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d_50 : float</span>
<span class="sd">        Volume median diameter (m)</span>
<span class="sd">    d_max : float</span>
<span class="sd">        Maximum stable diameter (m)</span>
<span class="sd">    alpha : float, default=1.8</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d_50 : float</span>
<span class="sd">        Volume median diameter (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function follows the idea of Sintef to not let d_95 of the Rosin-</span>
<span class="sd">    Rammler distribution exceed the maximum stable particle size.  If the </span>
<span class="sd">    original d_50 and d_max result in d_95 exceeding d_max, the d_50 is </span>
<span class="sd">    shifted downward such that d_95 will equal d_max.  Otherwise, the original</span>
<span class="sd">    d_50 is preserved.</span>
<span class="sd">    </span>
<span class="sd">    Uses the Rosin-Rammler distribution equations at:</span>
<span class="sd">    </span>
<span class="sd">        https://en.wikipedia.org/wiki/Particle-size_distribution</span>
<span class="sd">    </span>
<span class="sd">    last access on 03/17/20.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># k parameter for d50</span>
    <span class="n">k</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">0.5</span><span class="p">)</span>
    
    <span class="c1"># Adjust down if d95 exceeds the de_max</span>
    <span class="k">if</span> <span class="n">d_max</span> <span class="o">==</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="n">d50</span>
    
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Compute d95 for the given d50</span>
        <span class="n">d95</span> <span class="o">=</span> <span class="n">d50</span> <span class="o">*</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.95</span><span class="p">)</span> <span class="o">/</span> <span class="n">k</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">alpha</span><span class="p">)</span>
        
        <span class="c1"># Adjust d50 so that d95 does not exceed d_max</span>
        <span class="k">if</span> <span class="n">d95</span> <span class="o">&gt;</span> <span class="n">d_max</span><span class="p">:</span>
            <span class="n">d95</span> <span class="o">=</span> <span class="n">d_max</span>
            <span class="n">k95</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">0.05</span><span class="p">)</span>
            <span class="n">d50</span> <span class="o">=</span> <span class="n">d95</span> <span class="o">*</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="mf">0.5</span><span class="p">)</span> <span class="o">/</span> <span class="n">k95</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">alpha</span><span class="p">)</span>
    
    <span class="c1"># Return the final distribution fit</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span></div>


<div class="viewcode-block" id="log_normal_fit"><a class="viewcode-back" href="../autodoc/psd/psf.log_normal_fit.html#psf.log_normal_fit">[docs]</a><span class="k">def</span> <span class="nf">log_normal_fit</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">d_max</span><span class="p">,</span> <span class="n">sigma</span><span class="o">=</span><span class="mf">0.27</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return d_50 and sigma for the log-normal distribution</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d_50 : float</span>
<span class="sd">        Volume median diameter (m)</span>
<span class="sd">    d_max : float</span>
<span class="sd">        Maximum stable diameter (m)</span>
<span class="sd">    sigma : float, default=0.27</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d_50 : float</span>
<span class="sd">        Volume median diameter (m)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function follows the idea of Sintef to not let d_95 of the particle</span>
<span class="sd">    size distribution exceed the maximum stable particle size.  If the </span>
<span class="sd">    original d_50 and d_max result in d_95 exceeding d_max, the d_50 is </span>
<span class="sd">    shifted downward such that d_95 will equal d_max.  Otherwise, the original</span>
<span class="sd">    d_50 is preserved.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Adjust down if d50 exceeds the de_max</span>
    <span class="k">if</span> <span class="n">d_max</span> <span class="o">==</span> <span class="kc">None</span><span class="p">:</span>
        <span class="c1"># Do not adjust the fit</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="n">d50</span>
        
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Comnpute d95 for the given d50 and sigma</span>
        <span class="n">mu</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d50</span><span class="p">)</span>
        <span class="n">mu_95</span> <span class="o">=</span> <span class="n">mu</span> <span class="o">+</span> <span class="mf">1.6449</span> <span class="o">*</span> <span class="n">sigma</span>
        <span class="n">d95</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">mu_95</span><span class="p">)</span>
        
        <span class="c1"># Adjust d_50 so that d_95 does not exceed d_max</span>
        <span class="k">if</span> <span class="n">d95</span> <span class="o">&gt;</span> <span class="n">d_max</span><span class="p">:</span>
            <span class="n">d50</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">d_max</span><span class="p">)</span> <span class="o">-</span> <span class="mf">1.6449</span> <span class="o">*</span> <span class="n">sigma</span><span class="p">)</span>
    
    <span class="c1"># Return the final distribution fit</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span></div>


<span class="c1"># Functions for computing maximum stable particle size -----------------------</span>

<div class="viewcode-block" id="de_max_oil"><a class="viewcode-back" href="../autodoc/psd/psf.de_max_oil.html#psf.de_max_oil">[docs]</a><span class="k">def</span> <span class="nf">de_max_oil</span><span class="p">(</span><span class="n">rho_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculate the maximum stable oil droplet size</span>
<span class="sd">    </span>
<span class="sd">    Predicts the maximum stable liquid particle size per Clift et al. (1978) </span>
<span class="sd">    via the equation:</span>
<span class="sd">    </span>
<span class="sd">    d_max = 4. * np.sqrt(sigma / (g (rho - rho_p)))</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    rho_p : float</span>
<span class="sd">        Density of the phase undergoing breakup (kg/m^3)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Interfacial tension between the phase undergoing breakup and the </span>
<span class="sd">        ambient receiving continuous phase (N/m)</span>
<span class="sd">    rho : float</span>
<span class="sd">        Density of the ambient receiving continuous phase (kg/m^3) </span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d_max : float</span>
<span class="sd">        Maximum stable particle size (m)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">sigma</span> <span class="o">/</span> <span class="p">(</span><span class="n">G</span> <span class="o">*</span> <span class="p">(</span><span class="n">rho</span> <span class="o">-</span> <span class="n">rho_p</span><span class="p">)))</span></div>


<span class="k">def</span> <span class="nf">grow_rate</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">K</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the instability growth rate on a gas bubble</span>
<span class="sd">    </span>
<span class="sd">    Write instability growth rate equation in Grace et al. as a root</span>
<span class="sd">    problem for n = f(k)</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    res : float</span>
<span class="sd">        The residual of the growth-rate equation expressed as a root-finding</span>
<span class="sd">        problem.</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is used by the `grace()` function for maximum stable </span>
<span class="sd">    particle size.  It should not be called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute more derived variables</span>
    <span class="n">m_c</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">k</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">n</span> <span class="o">/</span> <span class="n">nu_c</span><span class="p">)</span>
    <span class="n">m_d</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">k</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">n</span> <span class="o">/</span> <span class="n">nu_d</span><span class="p">)</span>
    <span class="n">mu_c</span> <span class="o">=</span> <span class="n">nu_c</span> <span class="o">/</span> <span class="n">rho_c</span>
    
    <span class="c1"># Compute the residual of the root function</span>
    <span class="n">res</span> <span class="o">=</span> <span class="p">(</span><span class="n">sigma</span> <span class="o">*</span> <span class="n">k</span><span class="o">**</span><span class="mi">3</span> <span class="o">-</span> <span class="n">g</span> <span class="o">*</span> <span class="n">k</span> <span class="o">*</span> <span class="n">dp</span> <span class="o">+</span> <span class="n">n</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="p">(</span><span class="n">rho_c</span> <span class="o">+</span> <span class="n">rho_d</span><span class="p">))</span> <span class="o">*</span> \
          <span class="p">(</span><span class="n">k</span> <span class="o">+</span> <span class="n">m_c</span> <span class="o">+</span> <span class="n">K</span> <span class="o">*</span> <span class="p">(</span><span class="n">k</span> <span class="o">+</span> <span class="n">m_d</span><span class="p">))</span> <span class="o">+</span> <span class="mi">4</span> <span class="o">*</span> <span class="n">n</span> <span class="o">*</span> <span class="n">k</span> <span class="o">*</span> <span class="n">mu_c</span> <span class="o">*</span> <span class="p">(</span><span class="n">k</span> <span class="o">+</span> <span class="n">K</span> <span class="o">*</span> <span class="n">m_d</span><span class="p">)</span> <span class="o">*</span> \
          <span class="p">(</span><span class="n">K</span> <span class="o">*</span> <span class="n">k</span> <span class="o">+</span> <span class="n">m_c</span><span class="p">)</span>
    
    <span class="c1"># Return the residual</span>
    <span class="k">return</span> <span class="n">res</span>

<span class="k">def</span> <span class="nf">grow_time</span><span class="p">(</span><span class="n">lam</span><span class="p">,</span> <span class="n">de</span><span class="p">,</span> <span class="n">U</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">K</span><span class="p">,</span> <span class="n">c_0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compare the available and needed disturbance growth times for instability</span>
<span class="sd">    </span>
<span class="sd">    Compares the time available for a disturbance to grow to the time needed</span>
<span class="sd">    for that disturbance to break a bubble.  </span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    t_cr : float</span>
<span class="sd">        The critical time (s) for which the required grow time equals the</span>
<span class="sd">        available time</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is used by the `grace()` function for maximum stable </span>
<span class="sd">    particle size.  It should not be called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute the derived variables</span>
    <span class="n">k</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">/</span> <span class="n">lam</span>
    
    <span class="c1"># Consider disturbances with a node at the nose</span>
    <span class="n">theta_1</span> <span class="o">=</span> <span class="n">lam</span> <span class="o">/</span> <span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">de</span><span class="p">)</span>
    
    <span class="c1"># Compute the available time for disturbance growth</span>
    <span class="n">t_a</span> <span class="o">=</span> <span class="n">de</span> <span class="o">/</span> <span class="mi">2</span> <span class="o">/</span> <span class="n">U</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">+</span> <span class="mf">3.</span> <span class="o">/</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">K</span><span class="p">)</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="mf">1.</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">tan</span><span class="p">(</span>
          <span class="n">theta_1</span> <span class="o">/</span> <span class="mf">2.</span><span class="p">))</span>
    
    <span class="c1"># Compute the grwoth rate of this disturbance</span>
    <span class="n">n0</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="n">t_a</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">fsolve</span><span class="p">(</span><span class="n">grow_rate</span><span class="p">,</span> <span class="n">n0</span><span class="p">,</span> <span class="n">args</span><span class="o">=</span><span class="p">(</span><span class="n">k</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> 
                                 <span class="n">rho_d</span><span class="p">,</span> <span class="n">K</span><span class="p">)</span>
              <span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
    
    <span class="c1"># Relate n to t_e</span>
    <span class="n">t_e</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="n">n</span>
    
    <span class="c1"># Return the critical growth time</span>
    <span class="k">return</span> <span class="n">c_0</span> <span class="o">*</span> <span class="n">t_e</span> <span class="o">-</span> <span class="n">t_a</span>

<span class="k">def</span> <span class="nf">find_de</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">mu_d</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">K</span><span class="p">,</span> 
            <span class="n">lam_crit</span><span class="p">,</span> <span class="n">c_0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Search for the critical stable bubble size</span>
<span class="sd">    </span>
<span class="sd">    Search for the maximum stable bubble size of a gas bubble in water using</span>
<span class="sd">    the method in Grace et al.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    t_min : float</span>
<span class="sd">        The minimum time required for a disturbance of the given size to </span>
<span class="sd">        break the fluid particle.  </span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is used by the `grace()` function for maximum stable </span>
<span class="sd">    particle size.  It should not be called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Time available for growth, t_a</span>
    <span class="c1"># The travel time from the position where disturance starts to the equator</span>
    <span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">dbm_f</span>
    
    <span class="c1"># Compute the rise velocity of this bubble size</span>
    <span class="n">shape</span> <span class="o">=</span> <span class="n">dbm_f</span><span class="o">.</span><span class="n">particle_shape</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">shape</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">U</span> <span class="o">=</span> <span class="n">dbm_f</span><span class="o">.</span><span class="n">us_sphere</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">)</span>
    <span class="k">elif</span> <span class="n">shape</span> <span class="o">==</span> <span class="mi">2</span><span class="p">:</span>
        <span class="n">U</span> <span class="o">=</span> <span class="n">dbm_f</span><span class="o">.</span><span class="n">us_ellipsoid</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">mu_d</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">U</span> <span class="o">=</span> <span class="n">dbm_f</span><span class="o">.</span><span class="n">us_spherical_cap</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">)</span>
    
    <span class="c1"># lam_max is upper limit on leading interface disturbance size</span>
    <span class="n">lam_max</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">de</span> <span class="o">/</span> <span class="mf">2.</span>
    
    <span class="c1"># Find the wave length that corresponds to the maximum growth rate</span>
    <span class="n">lam</span> <span class="o">=</span> <span class="n">minimize</span><span class="p">(</span><span class="n">grow_time</span><span class="p">,</span> <span class="mf">1.001</span> <span class="o">*</span> <span class="n">lam_crit</span><span class="p">,</span> <span class="n">args</span><span class="o">=</span><span class="p">(</span><span class="n">de</span><span class="p">,</span> <span class="n">U</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span>
                                                      <span class="n">sigma</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> 
                                                      <span class="n">rho_d</span><span class="p">,</span> <span class="n">K</span><span class="p">,</span> <span class="n">c_0</span><span class="p">),</span>
                   <span class="n">bounds</span><span class="o">=</span><span class="p">[(</span><span class="n">lam_crit</span><span class="p">,</span> <span class="n">lam_max</span><span class="p">)]</span>
                   <span class="p">)</span><span class="o">.</span><span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">t_min</span> <span class="o">=</span> <span class="n">grow_time</span><span class="p">(</span><span class="n">lam</span><span class="p">,</span> <span class="n">de</span><span class="p">,</span> <span class="n">U</span><span class="p">,</span> <span class="n">nu_c</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span><span class="n">sigma</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">K</span><span class="p">,</span> 
                      <span class="n">c_0</span><span class="p">)</span>
    
    <span class="c1"># Return the growth time</span>
    <span class="k">return</span> <span class="n">t_min</span>


<div class="viewcode-block" id="grace"><a class="viewcode-back" href="../autodoc/psd/psf.grace.html#psf.grace">[docs]</a><span class="k">def</span> <span class="nf">grace</span><span class="p">(</span><span class="n">rho_c</span><span class="p">,</span> <span class="n">rho_d</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">,</span> <span class="n">mu_d</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">fp_type</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Implement the Grace et al. algorithm for maximum stable particle size</span>
<span class="sd">    </span>
<span class="sd">    Computes the maximum stable particle size of an immiscible particle </span>
<span class="sd">    rising in stagnant water following a method in Grace et al. </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    rho_c : float</span>
<span class="sd">        Density of the continuous-phase ambient fluid (kg/m^3)</span>
<span class="sd">    rho_d : float</span>
<span class="sd">        Density of the immiscible fluid particle subject to breakup (kg/m^3)</span>
<span class="sd">    mu_c : float</span>
<span class="sd">        Dynamic viscosity of the continuous-phase ambient fluid (Pa s)</span>
<span class="sd">    mu_d : float</span>
<span class="sd">        Dynamic viscosity of the immiscible fluid particle subject to breakup </span>
<span class="sd">        (Pa s)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Interfacial tension between the continuous phase ambient fluid and</span>
<span class="sd">        the immiscible fluid particle subject to breakup (N/m)</span>
<span class="sd">    fp_type : int, default=0</span>
<span class="sd">        Phase of the immiscible fluid particle; 0 = gas, 1 = liquid.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Equivalent spherical diameter of the maximum stable fluid particle</span>
<span class="sd">        subject to breakup in stagnant water (m)</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    grow_rate, grow_time, find_de</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    Implements the method in * Grace, J.R., Wairegi, T., Brophy, J., (1978)</span>
<span class="sd">    &quot;Break-up of drops and bubbles in stagnant media,&quot; Can. J. Chem. Eng. 56</span>
<span class="sd">    (1), 3-8.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    
    <span class="c1"># Set the fit parameter</span>
    <span class="k">if</span> <span class="n">fp_type</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="c1"># This is gas</span>
        <span class="n">c_0</span> <span class="o">=</span> <span class="mf">3.8</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># This is liquid</span>
        <span class="n">c_0</span> <span class="o">=</span> <span class="mf">1.4</span>
    
    <span class="c1"># Compute the derived properties</span>
    <span class="n">dp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">rho_c</span> <span class="o">-</span> <span class="n">rho_d</span><span class="p">)</span>
    <span class="n">K</span> <span class="o">=</span> <span class="n">mu_d</span> <span class="o">/</span> <span class="n">mu_c</span>
    <span class="n">nu_c</span> <span class="o">=</span> <span class="n">mu_c</span> <span class="o">/</span> <span class="n">rho_c</span>
    <span class="n">nu_d</span> <span class="o">=</span> <span class="n">mu_d</span> <span class="o">/</span> <span class="n">rho_d</span>
    
    <span class="c1"># Region of instability.</span>
    <span class="c1"># lam_crit is lower limit on unstable wavelengths</span>
    <span class="n">lam_crit</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">sigma</span> <span class="o">/</span> <span class="p">(</span><span class="n">G</span> <span class="o">*</span> <span class="n">dp</span><span class="p">))</span>
    
    <span class="c1"># Lower limit on maximum stable diameter is lam_crit = lam_max</span>
    <span class="n">de_max_star</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">lam_crit</span>
    
    <span class="c1"># Choose Initialize the search near this minimum</span>
    <span class="n">de</span> <span class="o">=</span> <span class="mf">1.01</span> <span class="o">*</span> <span class="n">de_max_star</span>
    
    <span class="c1"># Find the maximum stable bubble size</span>
    <span class="n">de_max</span> <span class="o">=</span> <span class="n">fsolve</span><span class="p">(</span><span class="n">find_de</span><span class="p">,</span> <span class="n">de</span><span class="p">,</span> <span class="n">args</span><span class="o">=</span><span class="p">(</span><span class="n">rho_d</span><span class="p">,</span> <span class="n">rho_c</span><span class="p">,</span> <span class="n">mu_d</span><span class="p">,</span> <span class="n">mu_c</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">nu_d</span><span class="p">,</span> 
                                       <span class="n">nu_c</span><span class="p">,</span> <span class="n">G</span><span class="p">,</span> <span class="n">dp</span><span class="p">,</span> <span class="n">K</span><span class="p">,</span> <span class="n">lam_crit</span><span class="p">,</span> <span class="n">c_0</span><span class="p">)</span>
                   <span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
    
    <span class="c1"># Return the result</span>
    <span class="k">return</span> <span class="n">de_max</span></div>


<span class="c1"># SINTEF Model Equations -----------------------------------------------------</span>

<div class="viewcode-block" id="sintef"><a class="viewcode-back" href="../autodoc/psd/psf.sintef.html#psf.sintef">[docs]</a><span class="k">def</span> <span class="nf">sintef</span><span class="p">(</span><span class="n">d0</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> 
           <span class="n">fp_type</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">use_d95</span><span class="o">=</span><span class="kc">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute characteristic values for jet breakup</span>
<span class="sd">    </span>
<span class="sd">    Computes the characteristic particle sizes for jet breakup using the</span>
<span class="sd">    equations in Johansen et al. (2013) (sintef model equations) and using</span>
<span class="sd">    the model coefficients updated in technical reports to API.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d0 : float</span>
<span class="sd">        Equivalent circular diameter of the release (m)</span>
<span class="sd">    m_gas : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the gas-phase fluid at the</span>
<span class="sd">        release (kg/s)</span>
<span class="sd">    rho_gas : float</span>
<span class="sd">        Density of the gas-phase petroleum fluid at the release (kg/m^3)</span>
<span class="sd">    m_oil : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the liquid-phase petroleum</span>
<span class="sd">        fluid at the release (kg/s)</span>
<span class="sd">    rho_oil : float</span>
<span class="sd">        Density of the liquid-phase petroleum fluid at the release (kg/m^3)</span>
<span class="sd">    mu_p : float</span>
<span class="sd">        Dynamic viscosity of the fluid phase of interest at the release </span>
<span class="sd">        (Pa s)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Interfacial tension between the fluid phase of interest and water at</span>
<span class="sd">        the release (N/m)</span>
<span class="sd">    rho : float</span>
<span class="sd">        Density of seawater at the release (kg/m^3)</span>
<span class="sd">    mu : float</span>
<span class="sd">        Dynamic viscosity of seawater at the release (Pa s)</span>
<span class="sd">    fp_type : int, default=1</span>
<span class="sd">        Fluid phase to compute breakup; 0 = gas, 1 = liquid.  The SINTEF</span>
<span class="sd">        equation authors do not recommend using this method for gas; hence,</span>
<span class="sd">        fp_type should normally equal 1.</span>
<span class="sd">    use_d95 : bool, default=True</span>
<span class="sd">        Flag indicating whether or not to implement the d_95 rule (see</span>
<span class="sd">        module documentation above); `True` means to use the rule.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Convert mass-flux to volume flux</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">m_gas</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">q_gas</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">q_gas</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">m_oil</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">q_oil</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">q_oil</span> <span class="o">=</span> <span class="mf">0.</span>
    
    <span class="c1"># Get the void-fraction adjusted characteristic exit velocity</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">q_gas</span> <span class="o">/</span> <span class="p">(</span><span class="n">q_gas</span> <span class="o">+</span> <span class="n">q_oil</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">q_oil</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
        <span class="c1"># This is gas only</span>
        <span class="n">Un</span> <span class="o">=</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">q_gas</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        <span class="n">rho_m</span> <span class="o">=</span> <span class="n">rho_gas</span>
    <span class="k">elif</span> <span class="n">q_gas</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="c1"># This is oil only</span>
        <span class="n">Un</span> <span class="o">=</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">q_oil</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        <span class="n">rho_m</span> <span class="o">=</span> <span class="n">rho_oil</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># This is oil and gas</span>
        <span class="n">Un</span> <span class="o">=</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">q_oil</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">n</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
        <span class="n">rho_m</span> <span class="o">=</span> <span class="n">rho_oil</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">n</span><span class="p">)</span> <span class="o">+</span> <span class="n">rho_gas</span> <span class="o">*</span> <span class="n">n</span>
    
    <span class="n">Fr</span> <span class="o">=</span> <span class="n">Un</span> <span class="o">/</span> <span class="p">(</span><span class="n">G</span> <span class="o">*</span> <span class="p">(</span><span class="n">rho</span> <span class="o">-</span> <span class="n">rho_m</span><span class="p">)</span> <span class="o">/</span> <span class="n">rho</span> <span class="o">*</span> <span class="n">d0</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
    <span class="n">Uc</span> <span class="o">=</span> <span class="n">Un</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">+</span> <span class="mf">1.</span><span class="o">/</span><span class="n">Fr</span><span class="p">)</span>
    
    <span class="c1"># Compute the particle size distribution parameters</span>
    <span class="k">if</span> <span class="n">fp_type</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">sintef_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> 
                                             <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> 
                                             <span class="n">use_d95</span><span class="o">=</span><span class="n">use_d95</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">sintef_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> 
                                             <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> 
                                             <span class="n">use_d95</span><span class="o">=</span><span class="n">use_d95</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span></div>


<span class="k">def</span> <span class="nf">sintef_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> 
                 <span class="n">use_d95</span><span class="o">=</span><span class="kc">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Computes the particle size for the Sintef equation</span>
<span class="sd">    </span>
<span class="sd">    Evaluates the parameters of the particle size distribution for the </span>
<span class="sd">    SINTEF equation and implements the d_95 rule as appropriate.  This</span>
<span class="sd">    function returns the parameters of the Rosin-Rammler size distribution</span>
<span class="sd">    with the spreading rates as reported in Johansen et al.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `sintef()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">if</span> <span class="n">q</span> <span class="o">&gt;</span> <span class="mf">0.</span><span class="p">:</span>
        
        <span class="c1"># Compute d_50 from the We model</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="n">sintef_d50</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">)</span>
        
        <span class="c1"># Get an estimate of de_max</span>
        <span class="k">if</span> <span class="n">is_gas</span><span class="p">:</span>
            <span class="n">de_max</span> <span class="o">=</span> <span class="n">grace</span><span class="p">(</span><span class="n">rho</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">fp_type</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">de_max</span> <span class="o">=</span> <span class="n">de_max_oil</span><span class="p">(</span><span class="n">rho_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">)</span>
        
        <span class="c1"># Get the adjusted particle size distribution</span>
        <span class="n">d50_from95</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">rosin_rammler_fit</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">)</span>
        
        <span class="c1"># Return the desired value for d50</span>
        <span class="k">if</span> <span class="n">use_d95</span><span class="p">:</span>
            <span class="c1"># Use the d_95 rule</span>
            <span class="n">d50</span> <span class="o">=</span> <span class="n">d50_from95</span>
        <span class="k">elif</span> <span class="n">d50</span> <span class="o">&gt;</span> <span class="n">de_max</span><span class="p">:</span>
            <span class="c1"># Truncate the distribution</span>
            <span class="n">d50</span> <span class="o">=</span> <span class="n">de_max</span>
    
    <span class="k">else</span><span class="p">:</span>
        
        <span class="c1"># Return an empty set of particles</span>
        <span class="n">de_max</span> <span class="o">=</span> <span class="kc">None</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">rosin_rammler_fit</span><span class="p">(</span><span class="mf">0.</span><span class="p">,</span> <span class="n">de_max</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span>


<span class="k">def</span> <span class="nf">sintef_d50</span><span class="p">(</span><span class="n">u0</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute d_50 from the SINTEF equations</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `sintef()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute the non-dimensional constants</span>
    <span class="n">We</span> <span class="o">=</span> <span class="n">rho_p</span> <span class="o">*</span> <span class="n">u0</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">d0</span> <span class="o">/</span> <span class="n">sigma</span>
    <span class="n">Vi</span> <span class="o">=</span> <span class="n">mu_p</span> <span class="o">*</span> <span class="n">u0</span> <span class="o">/</span> <span class="n">sigma</span>
    
    <span class="k">if</span> <span class="n">We</span> <span class="o">&gt;</span> <span class="mf">350.</span><span class="p">:</span>
        <span class="c1"># Atomization...use the the We model</span>
        <span class="n">A</span> <span class="o">=</span> <span class="mf">24.8</span>
        <span class="n">B</span> <span class="o">=</span> <span class="mf">0.08</span>
        
        <span class="c1"># Solve for the volume mean diameter from the implicit equation</span>
        <span class="k">def</span> <span class="nf">residual</span><span class="p">(</span><span class="n">dp</span><span class="p">):</span>
            <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">            Compute the residual of the SINTEF modified Weber number model</span>
<span class="sd">            </span>
<span class="sd">            Evaluate the residual of the non-dimensional diameter </span>
<span class="sd">            dp = de_50 / D for the SINTEF droplet break-up model.</span>
<span class="sd">            </span>
<span class="sd">            Input variables are:</span>
<span class="sd">                We, Vi, A, B = constant and global from above</span>
<span class="sd">                dp = Non-dimensional diameter de_50 / D (--)</span>
<span class="sd">            </span>
<span class="sd">            &quot;&quot;&quot;</span>
            <span class="c1"># Compute the non-dimensional diameter and return the residual</span>
            <span class="k">return</span> <span class="n">dp</span> <span class="o">-</span> <span class="n">A</span> <span class="o">*</span> <span class="p">(</span><span class="n">We</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">+</span> <span class="n">B</span> <span class="o">*</span> <span class="n">Vi</span> <span class="o">*</span> <span class="n">dp</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="mf">3.</span><span class="p">)))</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mf">3.</span><span class="o">/</span><span class="mf">5.</span><span class="p">)</span>
        
        <span class="c1"># Find the gas and liquid fraction for the mixture</span>
        <span class="n">dp</span> <span class="o">=</span> <span class="n">fsolve</span><span class="p">(</span><span class="n">residual</span><span class="p">,</span> <span class="mf">5.</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
        
        <span class="c1"># Compute the final d_50</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="n">dp</span> <span class="o">*</span> <span class="n">d0</span>
    
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Sinuous wave breakup...use the pipe diameter</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="mf">1.2</span> <span class="o">*</span> <span class="n">d0</span>
    
    <span class="c1"># Return the result</span>
    <span class="k">return</span> <span class="n">d50</span>


<span class="c1"># Li et al. Equations --------------------------------------------------------</span>

<div class="viewcode-block" id="li_etal"><a class="viewcode-back" href="../autodoc/psd/psf.li_etal.html#psf.li_etal">[docs]</a><span class="k">def</span> <span class="nf">li_etal</span><span class="p">(</span><span class="n">d0</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> 
            <span class="n">fp_type</span><span class="o">=</span><span class="mi">1</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute characteristic values for jet breakup</span>
<span class="sd">    </span>
<span class="sd">    Computes the characteristic particle sizes for jet breakup using the</span>
<span class="sd">    equations in Li et al. (2017) (li_etal model equations).  The authors</span>
<span class="sd">    provide different equation fit parameters for gas or liquid breakup, </span>
<span class="sd">    and this function selects the correct parameters for the fluid phase </span>
<span class="sd">    of interest.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d0 : float</span>
<span class="sd">        Equivalent circular diameter of the release (m)</span>
<span class="sd">    m_gas : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the gas-phase fluid at the</span>
<span class="sd">        release (kg/s)</span>
<span class="sd">    rho_gas : float</span>
<span class="sd">        Density of the gas-phase petroleum fluid at the release (kg/m^3)</span>
<span class="sd">    m_oil : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the liquid-phase petroleum</span>
<span class="sd">        fluid at the release (kg/s)</span>
<span class="sd">    rho_oil : float</span>
<span class="sd">        Density of the liquid-phase petroleum fluid at the release (kg/m^3)</span>
<span class="sd">    mu_p : float</span>
<span class="sd">        Dynamic viscosity of the fluid phase of interest at the release </span>
<span class="sd">        (Pa s)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Interfacial tension between the fluid phase of interest and water at</span>
<span class="sd">        the release (N/m)</span>
<span class="sd">    rho : float</span>
<span class="sd">        Density of seawater at the release (kg/m^3)</span>
<span class="sd">    mu : float</span>
<span class="sd">        Dynamic viscosity of seawater at the release (Pa s)</span>
<span class="sd">    fp_type : int, default=1</span>
<span class="sd">        Fluid phase to compute breakup; 0 = gas, 1 = liquid.  The SINTEF</span>
<span class="sd">        equation authors do not recommend using this method for gas; hence,</span>
<span class="sd">        fp_type should normally equal 1.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Convert mass-flux to volume flux</span>
    <span class="n">q_gas</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">)</span>
    <span class="n">q_oil</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">)</span>
    
    <span class="c1"># Get the void-fraction adjusted characteristic exit velocity</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">q_gas</span> <span class="o">/</span> <span class="p">(</span><span class="n">q_gas</span> <span class="o">+</span> <span class="n">q_oil</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">fp_type</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">Uc</span> <span class="o">=</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">q_gas</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">n</span>
    <span class="k">elif</span> <span class="n">fp_type</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">Uc</span> <span class="o">=</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">q_oil</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">n</span><span class="p">)</span>
    
    <span class="c1"># Compute the particle size distribution for gas and oil</span>
    <span class="k">if</span> <span class="n">fp_type</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">li_etal_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q_gas</span><span class="p">,</span> <span class="n">rho_gas</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> 
                                              <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">li_etal_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q_oil</span><span class="p">,</span> <span class="n">rho_oil</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> 
                                              <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span></div>


<span class="k">def</span> <span class="nf">li_etal_model</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Computes the particle size for the Li et al. equation</span>
<span class="sd">    </span>
<span class="sd">    Evaluates the parameters of the particle size distribution for the Li et</span>
<span class="sd">    al. equation and implements the d_95 rule. This function returns the</span>
<span class="sd">    parameters of the Rosin-Rammler size distribution with the spreading</span>
<span class="sd">    rates as reported in Li et al.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    k : float</span>
<span class="sd">        Scale parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    alpha : float</span>
<span class="sd">        Shape parameter for the Rosin-Rammler size distribution (--)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `li_etal()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">if</span> <span class="n">q</span> <span class="o">&gt;</span> <span class="mf">0.</span><span class="p">:</span>
        
        <span class="c1"># Compute d_50 from the We model</span>
        <span class="n">d50</span> <span class="o">=</span> <span class="n">li_etal_d50</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">is_gas</span><span class="p">)</span>
        
        <span class="c1"># Get an estimate of de_max</span>
        <span class="k">if</span> <span class="n">is_gas</span><span class="p">:</span>
            <span class="n">de_max</span> <span class="o">=</span> <span class="n">grace</span><span class="p">(</span><span class="n">rho</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">fp_type</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">de_max</span> <span class="o">=</span> <span class="n">de_max_oil</span><span class="p">(</span><span class="n">rho_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">)</span>
        
        <span class="c1"># Get the adjusted particle size distribution</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">rosin_rammler_fit</span><span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="kc">None</span><span class="p">)</span>
    
    <span class="k">else</span><span class="p">:</span>
        
        <span class="c1"># Return an empty set of particles</span>
        <span class="n">de_max</span> <span class="o">=</span> <span class="kc">None</span>
        <span class="n">d50</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">rosin_rammler_fit</span><span class="p">(</span><span class="mf">0.</span><span class="p">,</span> <span class="n">de_max</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">alpha</span><span class="p">)</span>


<span class="k">def</span> <span class="nf">li_etal_d50</span><span class="p">(</span><span class="n">Uc</span><span class="p">,</span> <span class="n">d0</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">mu_p</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">is_gas</span><span class="o">=</span><span class="kc">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute d_50 from the Li et al. equations</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the fluid phase of interest (m)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `li_etal()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the constants of the fit for a jet release</span>
    <span class="n">p</span> <span class="o">=</span> <span class="mf">0.460</span>
    <span class="n">q</span> <span class="o">=</span> <span class="o">-</span><span class="mf">0.518</span>
    <span class="k">if</span> <span class="n">is_gas</span><span class="p">:</span>
        <span class="n">r</span> <span class="o">=</span> <span class="mf">2.988</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">r</span> <span class="o">=</span> <span class="mf">14.05</span>
    
    <span class="c1"># The Li et al. (2017) model wrongly uses the oil equation for the </span>
    <span class="c1"># maximum stable particle size of gas</span>
    <span class="n">de_max</span> <span class="o">=</span> <span class="n">de_max_oil</span><span class="p">(</span><span class="n">sigma</span><span class="p">,</span> <span class="n">rho_p</span><span class="p">,</span> <span class="n">rho</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">de_max</span> <span class="o">&lt;</span> <span class="n">d0</span><span class="p">:</span>
        <span class="n">dc</span> <span class="o">=</span> <span class="n">de_max</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">dc</span> <span class="o">=</span> <span class="n">d0</span>
    
    <span class="c1"># Compute dimensionless groups</span>
    <span class="n">We</span> <span class="o">=</span> <span class="n">rho</span> <span class="o">*</span> <span class="n">Uc</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">dc</span> <span class="o">/</span> <span class="n">sigma</span>
    <span class="n">Oh</span> <span class="o">=</span> <span class="n">mu_p</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">rho_p</span> <span class="o">*</span> <span class="n">sigma</span> <span class="o">*</span> <span class="n">dc</span><span class="p">)</span>
    
    <span class="c1"># Compute the non-dimensional d50</span>
    <span class="n">ds</span> <span class="o">=</span> <span class="n">r</span> <span class="o">*</span> <span class="p">(</span><span class="mi">1</span> <span class="o">+</span> <span class="mi">10</span> <span class="o">*</span> <span class="n">Oh</span><span class="p">)</span><span class="o">**</span><span class="n">p</span> <span class="o">*</span> <span class="n">We</span><span class="o">**</span><span class="n">q</span>
    <span class="n">d50</span> <span class="o">=</span> <span class="n">ds</span> <span class="o">*</span> <span class="n">dc</span>
    
    <span class="k">return</span> <span class="n">d50</span>


<span class="c1"># Wang et al. Equations ------------------------------------------------------</span>

<div class="viewcode-block" id="wang_etal"><a class="viewcode-back" href="../autodoc/psd/psf.wang_etal.html#psf.wang_etal">[docs]</a><span class="k">def</span> <span class="nf">wang_etal</span><span class="p">(</span><span class="n">d0</span><span class="p">,</span> <span class="n">m_g</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> 
              <span class="n">m_l</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">rho_l</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">P</span><span class="o">=</span><span class="mf">4.e6</span><span class="p">,</span> <span class="n">T</span><span class="o">=</span><span class="mf">288.15</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute characteristic values for gas jet breakup</span>
<span class="sd">    </span>
<span class="sd">    Computes the characteristic gas bubble sizes for jet breakup using the</span>
<span class="sd">    equations in Wang et al. (2018) (wang_etal model equations).  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    d0 : float</span>
<span class="sd">        Equivalent circular diameter of the release (m)</span>
<span class="sd">    m_g : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the gas-phase fluid at the</span>
<span class="sd">        release (kg/s)</span>
<span class="sd">    rho_g : float</span>
<span class="sd">        Density of the gas-phase  fluid at the release (kg/m^3)</span>
<span class="sd">    mu_g : float</span>
<span class="sd">        Dynamic viscosity of the gas-phase fluid at the release (Pa s)</span>
<span class="sd">    sigma_g : float</span>
<span class="sd">        Interfacial tension between the gas-phase fluid and water at the</span>
<span class="sd">        release (N/m)</span>
<span class="sd">    rho : float</span>
<span class="sd">        Density of seawater at the release (kg/m^3)</span>
<span class="sd">    mu : float</span>
<span class="sd">        Dynamic viscosity of seawater at the release (Pa s)</span>
<span class="sd">    m_l : np.array</span>
<span class="sd">        Mass fluxes of each pseudo-component of the liquid-phase fluid at the </span>
<span class="sd">        release (kg/s)</span>
<span class="sd">    rho_l : float</span>
<span class="sd">        Density of the liquid-phase fluid at the release (kg/m^3)</span>
<span class="sd">    P : float, default=4.e6</span>
<span class="sd">        Pressure in the receiving fluid (Pa); used to compute the speed of</span>
<span class="sd">        sound in the released gas.</span>
<span class="sd">    T : float, default=288.15</span>
<span class="sd">        Temperature of the gas phase at the release (K); used to compute the</span>
<span class="sd">        speed of sound in the released gas.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50_gas : float</span>
<span class="sd">        Volume median diameter of the gas bubbles (m)</span>
<span class="sd">    m_gas : float</span>
<span class="sd">        Mass fluxes of each pseudo-component of the gas-phase fluid at the</span>
<span class="sd">        release (kg/s).  This may be different from the input value in the </span>
<span class="sd">        case of choked flow at the orifice.</span>
<span class="sd">    m_oil : float</span>
<span class="sd">        Mass fluxes of each pseudo-component of the liquid-phase fluid at the</span>
<span class="sd">        release (kg/s).  This may be different from the input value in the </span>
<span class="sd">        case of choked flow at the orifice.</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    sigma : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Convert mass-flux to volume flux</span>
    <span class="n">Qg</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_g</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">m_l</span><span class="p">)</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
        <span class="n">Ql</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">Ql</span> <span class="o">=</span> <span class="n">mass2vol</span><span class="p">(</span><span class="n">m_l</span><span class="p">,</span> <span class="n">rho_l</span><span class="p">)</span>
    
    <span class="c1"># Compute the exit velocity assuming no choked flow and single exit</span>
    <span class="c1"># velocity</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">Qg</span> <span class="o">/</span> <span class="p">(</span><span class="n">Qg</span> <span class="o">+</span> <span class="n">Ql</span><span class="p">)</span>
    <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">d0</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mf">4.</span>
    <span class="n">Ug</span> <span class="o">=</span> <span class="p">(</span><span class="n">Qg</span> <span class="o">+</span> <span class="n">Ql</span><span class="p">)</span> <span class="o">/</span> <span class="n">A</span>
    
    <span class="c1"># Check for choked flow using methane for speed of sound</span>
    <span class="n">ch4</span> <span class="o">=</span> <span class="n">dbm</span><span class="o">.</span><span class="n">FluidMixture</span><span class="p">([</span><span class="s1">&#39;methane&#39;</span><span class="p">])</span>
    <span class="n">delta_rho</span> <span class="o">=</span> <span class="n">ch4</span><span class="o">.</span><span class="n">density</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">1.</span><span class="p">]),</span> <span class="n">T</span><span class="p">,</span> <span class="n">P</span><span class="p">)[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> \
                <span class="n">ch4</span><span class="o">.</span><span class="n">density</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">1.</span><span class="p">]),</span> <span class="n">T</span><span class="p">,</span> <span class="mf">1.01</span> <span class="o">*</span> <span class="n">P</span><span class="p">)[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">P</span> <span class="o">-</span> <span class="mf">1.01</span> <span class="o">*</span> <span class="n">P</span><span class="p">)</span> <span class="o">/</span> <span class="n">delta_rho</span><span class="p">)</span>
    <span class="k">if</span> <span class="mf">10.</span> <span class="o">*</span> <span class="n">Ug</span> <span class="o">&lt;</span> <span class="n">a</span><span class="p">:</span>
        <span class="n">U_E</span> <span class="o">=</span> <span class="n">Ug</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Compute the cp / cv ratio</span>
        <span class="n">cp_ch4</span> <span class="o">=</span> <span class="mf">35.69</span>  <span class="c1"># J/mol/K;  CO2 = 37.13</span>
        <span class="n">cv_ch4</span> <span class="o">=</span> <span class="n">cp_ch4</span> <span class="o">-</span> <span class="mf">8.31451</span>  <span class="c1"># From Poling et al. for ideal gases</span>
        <span class="n">kappa</span> <span class="o">=</span> <span class="n">cp_ch4</span> <span class="o">/</span> <span class="n">cv_ch4</span>    <span class="c1"># Assume approximately ok for petroleum</span>
        
        <span class="c1"># Get the Mach number</span>
        <span class="n">Ma</span> <span class="o">=</span> <span class="n">Ug</span> <span class="o">/</span> <span class="n">a</span>
        
        <span class="c1"># Correct the exit velocity for choked flow</span>
        <span class="k">if</span> <span class="n">Ma</span> <span class="o">&lt;</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">kappa</span> <span class="o">+</span> <span class="mf">1.</span><span class="p">)</span> <span class="o">/</span> <span class="mf">2.</span><span class="p">):</span>
            <span class="n">U_E</span> <span class="o">=</span> <span class="n">a</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="mf">1.</span> <span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">1.</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="p">(</span><span class="n">kappa</span> <span class="o">-</span> <span class="mf">1.</span><span class="p">)</span> <span class="o">*</span> <span class="n">Ma</span><span class="o">**</span><span class="mf">2.</span><span class="p">))</span> <span class="o">/</span> \
                  <span class="p">((</span><span class="n">kappa</span> <span class="o">-</span> <span class="mf">1.</span><span class="p">)</span> <span class="o">*</span> <span class="n">Ma</span><span class="p">)</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">U_E</span> <span class="o">=</span> <span class="n">a</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span> <span class="o">/</span> <span class="p">(</span><span class="n">kappa</span> <span class="o">+</span> <span class="mf">1.</span><span class="p">))</span>
    
    <span class="c1"># Update the gas and oil exit velocities</span>
    <span class="k">if</span> <span class="n">Qg</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">Ug</span> <span class="o">=</span> <span class="n">U_E</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">Ug</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">if</span> <span class="n">Ql</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">Ul</span> <span class="o">=</span> <span class="n">U_E</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">Ul</span> <span class="o">=</span> <span class="mi">0</span>
    
    <span class="c1"># Compute the particle size distribution for gas</span>
    <span class="n">d50_gas</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">sigma</span> <span class="o">=</span> <span class="n">wang_etal_model</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">Ug</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> 
                                                           <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">Ul</span><span class="p">,</span> 
                                                           <span class="n">rho_l</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50_gas</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span></div>


<span class="k">def</span> <span class="nf">wang_etal_model</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">Ug</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">Ul</span><span class="p">,</span> <span class="n">rho_l</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Computes the particle size for the Wang et al. equation</span>
<span class="sd">    </span>
<span class="sd">    Evaluates the parameters of the gas bubble size distribution for the Wang</span>
<span class="sd">    et al. equation and implements the d_95 rule. This function returns the</span>
<span class="sd">    parameters of the log-normal size distribution with the spreading</span>
<span class="sd">    rates as reported in Wang et al.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50_gas : float</span>
<span class="sd">        Volume median diameter of the gas bubbles (m)</span>
<span class="sd">    m_gas : float</span>
<span class="sd">        Mass fluxes of each pseudo-component of the gas-phase fluid at the</span>
<span class="sd">        release (kg/s).  This may be different from the input value in the </span>
<span class="sd">        case of choked flow at the orifice.</span>
<span class="sd">    m_oil : float</span>
<span class="sd">        Mass fluxes of each pseudo-component of the liquid-phase fluid at the</span>
<span class="sd">        release (kg/s).  This may be different from the input value in the </span>
<span class="sd">        case of choked flow at the orifice.</span>
<span class="sd">    de_max : float</span>
<span class="sd">        Maximum stable particle size of the fluid phase of interest (m)</span>
<span class="sd">    sigma_ln : float</span>
<span class="sd">        Standard deviation of the Log-normal distribution in logarithmic </span>
<span class="sd">        units.</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `wang_etal()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">if</span> <span class="n">Ug</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        
        <span class="c1"># Compute d50 from the model</span>
        <span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">)</span> <span class="o">=</span> <span class="n">wang_etal_d50</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">Ug</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">Ul</span><span class="p">,</span> 
                                      <span class="n">rho_l</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">)</span>
        
        <span class="c1"># Compute the maximum stable bubble size</span>
        <span class="n">de_max</span> <span class="o">=</span> <span class="n">grace</span><span class="p">(</span><span class="n">rho</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> <span class="n">mu</span><span class="p">,</span> <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">fp_type</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
        
        <span class="c1"># Get the adjusted particle size distribution</span>
        <span class="n">d50_gas</span><span class="p">,</span> <span class="n">sigma_ln</span> <span class="o">=</span> <span class="n">log_normal_fit</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">sigma</span><span class="o">=</span><span class="mf">0.27</span><span class="p">)</span>
        
    <span class="k">else</span><span class="p">:</span>
        
        <span class="c1"># Return an empty set of particles</span>
        <span class="n">m_gas</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">m_oil</span> <span class="o">=</span> <span class="n">rho_l</span> <span class="o">*</span> <span class="n">A</span> <span class="o">*</span> <span class="n">Ul</span>
        <span class="n">de_max</span> <span class="o">=</span> <span class="kc">None</span>
        <span class="n">d50_gas</span><span class="p">,</span> <span class="n">sigma_ln</span> <span class="o">=</span> <span class="n">log_normal_fit</span><span class="p">(</span><span class="mf">0.</span><span class="p">,</span> <span class="n">de_max</span><span class="p">)</span>
    
    <span class="k">return</span> <span class="p">(</span><span class="n">d50_gas</span><span class="p">,</span> <span class="n">m_gas</span><span class="p">,</span> <span class="n">m_oil</span><span class="p">,</span> <span class="n">de_max</span><span class="p">,</span> <span class="n">sigma_ln</span><span class="p">)</span>


<span class="k">def</span> <span class="nf">wang_etal_d50</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">Ug</span><span class="p">,</span> <span class="n">rho_g</span><span class="p">,</span> <span class="n">mu_g</span><span class="p">,</span> <span class="n">sigma_g</span><span class="p">,</span> <span class="n">Ul</span><span class="p">,</span> <span class="n">rho_l</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">mu</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute d_50 from the Wang et al. equations</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    d50 : float</span>
<span class="sd">        Volume median diameter of the gas bubbles (m)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    This function is called by the `wang_etal()` function after several </span>
<span class="sd">    intermediate parameters are computed.  This function should not be </span>
<span class="sd">    called directly.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute the total dynamic momentum and buoyancy fluxes</span>
    <span class="n">Ag</span> <span class="o">=</span> <span class="n">A</span> <span class="o">*</span> <span class="n">n</span>
    <span class="n">Al</span> <span class="o">=</span> <span class="n">A</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">n</span><span class="p">)</span>
    <span class="n">mg</span> <span class="o">=</span> <span class="n">rho_g</span> <span class="o">*</span> <span class="n">Ag</span> <span class="o">*</span> <span class="n">Ug</span><span class="o">**</span><span class="mi">2</span>
    <span class="n">bg</span> <span class="o">=</span> <span class="p">(</span><span class="n">rho</span> <span class="o">-</span> <span class="n">rho_g</span><span class="p">)</span> <span class="o">*</span> <span class="n">G</span> <span class="o">*</span> <span class="n">Ag</span> <span class="o">*</span> <span class="n">Ug</span>
    <span class="k">if</span> <span class="n">n</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">ml</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">bl</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">ml</span> <span class="o">=</span> <span class="n">rho_l</span> <span class="o">*</span> <span class="n">Al</span> <span class="o">*</span> <span class="n">Ul</span><span class="o">**</span><span class="mi">2</span>
        <span class="n">bl</span> <span class="o">=</span> <span class="p">(</span><span class="n">rho</span> <span class="o">-</span> <span class="n">rho_l</span><span class="p">)</span> <span class="o">*</span> <span class="n">G</span> <span class="o">*</span> <span class="n">Al</span> <span class="o">*</span> <span class="n">Ul</span>
    <span class="n">mo</span> <span class="o">=</span> <span class="n">mg</span> <span class="o">+</span> <span class="n">ml</span>
    <span class="n">bo</span> <span class="o">=</span> <span class="n">bg</span> <span class="o">+</span> <span class="n">bl</span>
    
    <span class="c1"># The kinematic momentum and buoyancy fluxes are</span>
    <span class="n">M</span> <span class="o">=</span> <span class="n">mo</span> <span class="o">/</span> <span class="n">rho</span>
    <span class="n">B</span> <span class="o">=</span> <span class="n">bo</span> <span class="o">/</span> <span class="n">rho</span>
    
    <span class="c1"># Jet-to-plume transition length scale</span>
    <span class="n">l_M</span> <span class="o">=</span> <span class="n">M</span><span class="o">**</span><span class="p">(</span><span class="mf">3.</span><span class="o">/</span><span class="mf">4.</span><span class="p">)</span> <span class="o">/</span> <span class="n">B</span><span class="o">**</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
    
    <span class="c1"># Characteristic velocity scale</span>
    <span class="n">Ua</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">mo</span> <span class="o">/</span> <span class="p">(</span><span class="n">rho</span> <span class="o">*</span> <span class="n">A</span><span class="p">))</span>
    
    <span class="c1"># Compute the mixture density</span>
    <span class="k">if</span> <span class="n">n</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">rho_l</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="n">rho_m</span> <span class="o">=</span> <span class="n">n</span> <span class="o">*</span> <span class="n">rho_g</span> <span class="o">+</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">n</span><span class="p">)</span> <span class="o">*</span> <span class="n">rho_l</span>
    
    <span class="c1"># Get the modified Weber number</span>
    <span class="n">We_m</span> <span class="o">=</span> <span class="n">rho_m</span> <span class="o">*</span> <span class="n">Ua</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">l_M</span> <span class="o">/</span> <span class="n">sigma_g</span>
    
    <span class="c1"># Compute the characteristic droplet size</span>
    <span class="n">d</span> <span class="o">=</span> <span class="mf">4.3</span> <span class="o">*</span> <span class="n">We_m</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mf">3.</span><span class="o">/</span><span class="mf">5.</span><span class="p">)</span> <span class="o">*</span> <span class="n">l_M</span>
        
    <span class="c1"># Compute the actual gas and oil flow rate</span>
    <span class="n">m_g</span> <span class="o">=</span> <span class="n">rho_g</span> <span class="o">*</span> <span class="n">Ag</span> <span class="o">*</span> <span class="n">Ug</span>
    <span class="n">m_l</span> <span class="o">=</span> <span class="n">rho_l</span> <span class="o">*</span> <span class="n">Al</span> <span class="o">*</span> <span class="n">Ul</span>
    
    <span class="c1"># Return the characteristic size</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">m_g</span><span class="p">,</span> <span class="n">m_l</span><span class="p">)</span>
</pre></div>

          </div>
          
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
<h1 class="logo"><a href="../index.html">Texas A&M Oil spill / Outfall Calculator</a></h1>








<h3>Navigation</h3>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../user_manual.html">TAMOC User Manual</a></li>
<li class="toctree-l1"><a class="reference internal" href="../unit_tests.html">Unit Tests</a></li>
<li class="toctree-l1"><a class="reference internal" href="../glossary.html">Glossary</a></li>
</ul>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../bugs.html">Reporting Bugs</a></li>
<li class="toctree-l1"><a class="reference internal" href="../readme.html">Read Me File</a></li>
<li class="toctree-l1"><a class="reference internal" href="../release.html">Release Notes</a></li>
<li class="toctree-l1"><a class="reference internal" href="../license.html">License</a></li>
</ul>

<div class="relations">
<h3>Related Topics</h3>
<ul>
  <li><a href="../index.html">Documentation overview</a><ul>
  <li><a href="index.html">Module code</a><ul>
  </ul></li>
  </ul></li>
</ul>
</div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="../search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" />
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>$('#searchbox').show(0);</script>








        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="footer">
      &copy;2020, Scott A. Socolofsky.
      
      |
      Powered by <a href="http://sphinx-doc.org/">Sphinx 2.4.4</a>
      &amp; <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
      
    </div>

    

    
  </body>
</html>